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The scaling of the conductivity at the superfiuid-insulator quantum phase transition in two di- 
mensions is studied by numerical simulations of the Bose-Hubbard model. In contrast to previous 
studies, we focus on properties of this model in the experimentally relevant thermodynamic limit 
at finite temperature T. We find clear evidence for deviations from ajfc-scaling of the conductivity 
towards cuk /T-scaling at low Matsubara frequencies us^. By careful analytic continuation using Pade 
approximants we show that this behavior carries over to the real frequency axis where the conduc- 
tivity scales with us/T at small frequencies and low temperatures. We estimate the universal dc 
conductivity to be a* — 0.45(5)Q 2 fh, distinct from previous estimates in the T = 0, us/T ^> 1 limit. 
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The non-trivial properties of materials in the vicinity 
of quantum phase transitions (QPTs) are an object of 
intense theoretical and experimental studies. The 

effect of quantum fluctuations driving the QPTs is espe- 
cially pronounced in low-dimensional systems, such as 
high-temperature superconductors and two-dimensional 
(2D) electron gases, exhibiting the quantum Hall effect. 
Particularly valuable are theoretical predictions of the 
behavior of the dynamical response functions, such as 
the optical conductivity and the dynamic structure fac- 
tor, since they allow for direct comparison of the theoreti- 
cal results with experimental data. It was pointed out by 
Damle and Sachdev _2| , that at the quantum-critical cou- 
pling the scaled dynamic conductivity T^ 2 ~ d ^ z a(us, T) at 
low frequencies and temperatures is a function of the sin- 
gle variable Hus/ksT: 

cr(w/T, T ^ 0) = (k B T/hc) {d ~ 2)/z a Q H(Twj/k B T). (1) 

Here <tq = Q 2 /h is the conductivity "quantum" (Q = 2e 
for the models we consider), S(x = hus/ksT) is a uni- 
versal dimensionless scaling function, c a non-universal 
constant, and z the dynamical critical exponent. For 
d = 2 the exponent vanishes, leading to a purely uni- 
versal conductivity 0, depending only on frequency us, 
measured against a characteristic time h/3 set by finite 
temperature T as tuj/ksT. Once hus/ksT 3> 1, for fixed 
T, the system no longer "feels" the effect of finite tem- 
perature and it is natural to expect that at such high 
us a crossover to a temperature-independent regime will 
take place ||, so that a(us,T) ~ cr(u;) with cr(uj) de- 
caying at high frequencies as 1/us 2 0- Deviations from 
scaling of a with us therefore signal that temperature ef- 
fects have become important. Note that the predicted 
universal behavior occurs for fixed us/T as T — » 0. The 
physical mechanisms of transport are predicted 0] to be 
quite distinct in the different regimes determined by the 
value of the scaling variable x: hydrodynamic, collision- 
dominated for i C 1, and collisionlcss, phase-coherent 
for x 3> 1 with £ = S(oo) largely independent of x in 
d = 2 and a independent of T |2j, |£j . 
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FIG. 1: Mean-field ground state phase diagram of the 2D 
Bose-Hubbard model. Shaded areas are the Mott insulating 
phases. Dashed lines are constant density profiles in steps of 
0.2. • indicates the location of QCP at the tip of Mott lobe 
as determined by SSE simulations along the dotted line. 

Intriguingly, early numerical studies H H 13 of 
QPTs in model systems have failed to observe scaling 
with hus/ksT. The results of the experiments seeking 
to verify the scaling hypothesis are ambiguous as well. 
Some of them, performed at the 2D quantum Hall transi- 
tions 0] an d 3D metal- insulator transitions , appear 
to s upp ort it. Others either note the absence of scal- 
ing 01 or suggest a different scaling form 0]. While 
the discrepancy between theory and experiment may be 
attributed to the unsuitable choice of the measurement 
regime 2], typically leading to fiLo/ksT ^> 1, there is 
no good reason why the predicted scaling would not be 
observable in numerical simulations if careful extrapola- 
tions first to L — > oo and then T — > for fixed us/T are 
performed. 

Our primary goal is to resolve this controversy by per- 
forming precise numerical simulations of the frequency- 
dependent conductivity at finite temperatures in the 
vicinity of the 2D QPT, exploiting recent algorithmic ad- 
vances to access larger system sizes and a wider temper- 
ature range. After the extrapolation of the results to the 
thermodynamic and T = limits and careful analytic 
continuation, we are able to demonstrate how the pre- 
dicted universal behavior of the conductivity may indeed 
be revealed. 
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We consider the 2D Bose-Hubbard model with the 
Hamiltonian TCbh = T~to + 7~ti, where the first term de- 
scribes the non-interacting softcore bosons hopping via 
the nearest-neighbor links of a 2D square lattice, and the 
second one includes the Hubbard-like on-site interactions: 

n o = + bl +s K) - A*5>r, (2) 



n x = |$> r (n r -l). (3) 



Here 6 = x, y, n r — b\b T is the particle number operator 
on site r, and b\,b r are the boson creation and annihi- 
lation operators at site r. Model parameters include the 
hopping constant t, Hubbard repulsion U , and chemi- 
cal potential p. The mean-field ground state phase dia- 
gram of this model (Fig. displays a number of Mott- 
insulating lobes with fixed integer boson density at low 
Zt/U (Z = 4 is the lattice coordination number). As 
the hopping t is increased or p is varied, a QPT to a 
superfluid (SF) phase takes place. We concentrate on 
the QPT occuring at the tip of the Mott lobe along the 
path of constant p, distinct from the generic transition 
occurring elsewhere along the phase boundary |l4j ] . 

The numerical simulations of Hbh were performed us- 
ing the stochastic series expansion (SSE) technique with 
directed loop updates 0, 0] , which is known to be very 
efficient for the simulations of boson models. Further- 
more, as described below, it allows us to directly evaluate 
the relevant correlation functions without discretization 
or numerical integration over the imaginary time. We 
have also employed an alternative (2 + l)-dimensional 
classical representation 0, of Hbh in terms of link- 
current variables describing the total bosonic current 
J = (J x , J y , J T ) defined on a discrete L x L x L T space- 
time lattice (L T A T = ft/3): 



K ^ 



2 J (r,r) 



(4) 



J has to be conserved and is therefore divergence-free, 
V • J = 0. The link-current variables take on integer 
values J x ' y ' T = 0,±1,±2, ... and denote the deviation 
of the particle number from its mean, so the transition 
corresponds to p = 0. K is the effective temperature, 
varying like t/U in Hbh- The model defined by TLv has 
been studied in the past using a very efficient directed 
geometrical worm algorithm Jl8| , and its critical point at 
p = has been determined 18] to be K c = 0.33305(5). 
A drawback of this representation is that the time direc- 
tion is discrete, imposing an ultra-violet frequency cut-off 
of data at ui c — 1/Ar, in contrast to SSE, where there 
is no such problem. The two numerical approaches are 
therefore largely complementary. These advanced tech- 
niques allowed us to simulate lattices of linear sizes up 
to L = 30 and inverse temperatures up to (3 = 10 using 
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FIG. 2: The compressibility (a) and SF density (b) versus 
Zt/U at fixed fi/U = 0.375 and aspect ratio f3L~ z = 0.5. 
The dashed vertical line is drawn at Zt c /U = 0.2385. Error 
bars are displayed only if larger than the symbol size. 



SSE for Hbh, and L = 256, L T = 64 using the directed 
geometrical worm algorithm for Tiv ■ 

Performing SSE simulations of Hbb we first precisely 
locate the quantum critical point (QCP) for fixed fi/U 
at the tip of the first Mott lobe. This transition belongs 
to the (2+l)D XY universality class with a dynamical 
critical exponent z = 1 14] . In the vicinity of the QCP 
the SF density p s and compressibility k are expected to 
obey the scaling relations 



Ps = L 2 ~ d ~ z Yx{5L L/l ' , f3L~ z ), 
k = L z - d Y 2 (5L l l v ,(3L~ Z ). 



-l/v 



(•5) 
(6) 



Here v is the correlation length critical exponent, [3 = 
1/ksT, S = \t — t c \, and Yi^ix, y) are the two-variable 
scaling functions. For a fixed aspect ratio (3L~ Z plots 
of the Lp a and Lk should then intersect at the critical 
point 6 = 0. Results of such a calculation at constant 
p c /U — 0.375 are presented in Fig. |2 from which we 
determine Zt c /U = 0.2385(5), v = 0.66(5). The position 
of the QCP and the shape of the phase boundary in its 
vicinity is consistent with previous simulations [l9( and 
strong-coupling perturbation theory [2(| • 

To analyze the behavior of the zero-momentum con- 
ductivity in the vicinity of the QCP we employ the re- 
lation between the dynamic conductivity cr(u>) and the 
Fourier transform A xx (lo) of the time-dependent current- 
current correlation function (CCCF) established by the 
Kubo formula (HH^. In SSE the real-time CCCF re- 
quired to determine A xx (tu) is not directly accessible. In- 
stead, the standard approach is to measure the CCCF 
A xx (t) — (jx(T)j x (0)} on the imaginary time axis, cal- 
culate its Fourier transform A xx (iujk) as a function of 
the Matsubara frequencies lo^ = 2nk //3, and analytically 
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continue the result to real frequencies |a,|2i!22|- Here and 
below we adopt a unit system in which both Q and h are 
unity, and j x (r) is the Heisenberg representation of the 
current operator j x = it[bl +x b r — blb r+x \. We have: 



i ■ \ (-M - k xx (iu k ) _ p(im k ) 
<T(iuJ k ) = 2it<jq = 2-koq . (7) 



L0 k 



L0 k 



Here (—k x ) is the kinetic energy per link and p{iLo k ) is 
the frequency-dependent stiffness. To measure A xx (iu k ) 
we note that A xx (t) may be expressed in terms of the 
correlation functions A^(r,r) = (K2(r, t)K%(0, 0)) of 
operators K+(t,t) — t for +x (r)6 r (r) and K x (r,r) = 
1 6J(r)&r+ x (r), which may be estimated efficiently in 
SSE |l5(. Remarkably it is possible to analytically per- 
form the Fourier transform with respect to r yielding 




(u> k )N(v^\m) 



(8) 



where N(v, 7; m) is the number of times the operators 
iT 7 (r) and K u (0) appear in the SSE operator sequence 
separated by m operator positions, and n is the expansion 
order. The coefficients a mn {uj k ) are given by the degen- 
erate hypergeometric (Kummer) function: a mn (uj k ) — 
\Fi(m + l,n; —i(3uj k ). This expression and (jHJ allow us 
to directly evaluate A xx (r,u> k ) as a function of Matsub- 
ara frequencies, eliminating any errors associated with 
the discretization of the imaginary time interval. Anal- 
ogously, in the link-current representation p(iuj k ) can be 
calculated , and the conductivity can be obtained from 
Eq. Q. 

In Fig. 13 we show results for a{iuj k ) versus uj k ob- 
tained using the geometrical worm algorithm on TLy at 
K c (Fig. and by SSE simulations at t c , p c of Hbh 
(Fig. Efc) . In both cases the results have been extrapo- 
lated to the thermodynamic limit L — > 00 at fixed (3. As 
evident from Fig.|^, the results deviate from scaling with 
uj k at small uj k and more significantly so at higher tem- 
peratures (small L T ). These deviations are also visible in 
the continuous time SSE data in Fig. 0;, demonstrating 
that they cannot be attributed to time discretization er- 
rors. Similar deviations have been noted previously 0,0 
but were not analyzed at fixed f3. Since the deviations 
persist in the L — > 00 limit at fixed (3, they may only 
be interpreted as finite T effects. Expecting a crossover 
to u> k /T scaling at small tu k , we plot our results versus 
Lo k /T in Fig. 01. For L T > 32, <t(ll>i/T) is already inde- 
pendent of T (L T ). In fact, as shown in Fig.|2tl, for W1...5, 
a(u> k /T, T) can unambigously be extrapolated to a finite 
a(oj k /T,T — > 0) ~ T,(x) limit. This fact is a clear indi- 
cation that cjfe/T-scaling indeed occurs as T — » 0. Ten- 
tatively for increasing w k /T, a(uj k /T,T — > 0) appears 
to reach a constant value of roughly 0.33ctq ~ E(oo) 
in excellent agreement with theoretical estimates . 
We note that deviations from w^-scaling appear to be 



largely absent in simulations of Hbh with disorder 00. 
However, at this QCP the dynamical critical exponent 
is different (z — 2). As is evident from the size of the 
error bars in Fig. |3| simulations of Tty are much more 
efficient than the SSE simulations directly on 7Ybh- In 
the following analytic continuation we therefore use the 
SSE data mostly as a consistency check. 

Our results on the imaginary frequency axis are lim- 
ited by the lowest Matsubara frequency u>i — inksT/h. 
However, the information about the behavior of t'(uj) = 
Recr(w) at low u> is embedded in values of the CCCF at 
all Matsubara frequencies, allowing us to determine it. In 
order to study the w/T-scaling predicted for the hydro- 
dynamic collision dominated regime Cl,we 
have attempted analytic continuations of p(iuj k ) to obtain 
o~'(uj) at real frequencies. SSE results for TIbh were an- 
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FIG. 3: The conductivity a(u> k ) in units of <tq versus Matsub- 
ara frequency ui k /uj c as obtained from Tiv (a). All results have 
been extrapolated to the thermodynamic limit L — * 00 using 
the scaling form f(L) — a + 6exp(— L/Q/y/L [23| by calculat- 
ing p(ufc) at fixed L T using 9 lattice sizes from L — L T . . . 4L T 
as shown in (b). a(ui k ) in units of <tq versus Matsubara fre- 
quency u) k as obtained from SSE calculations of Hbh, with 
some typical error bars shown. All results have been extrapo- 
lated to the thermodynamic limit by calculating A xx (u>k) for 
fixed P using 5 lattice sizes L = 12 . . . 30 (c). Scaling plot 
of the conductivity data from (a) versus ui k /T. • denotes 
extrapolations to T — ► (L T —> 00) at fixed u) k /T using: 
f{L T ) = c + dexp(-L T /£ T )/v^7 23| (d). 
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FIG. 4: The real part of the conductivity a' at the critical coupling in units of <jq. The data marked L T , plotted versus uj/lj c , 
were obtained using Tiv , combined with the analytic continuation of p(u/u) c ) as explained in the text. Results for ui/ui c ^ 1/2 
are denoted by dotted lines. The data marked SSE, plotted versus w/10, were obtained by direct SSE simulations of TCbh with 
L = 20, /3 — 10 and subsequent maximum entropy analysis (a). Results as a function of ui/T (b). 



alytically continued using the Bryan maximum entropy 
(ME) method |25j with flat initial image. For the results 
obtained for the link-current model Tly we use a method 
that should be most sensitive to low frequencies lu/lo c < 1 
or T,(x <C 1). We fit the extrapolated low frequency part 
(first 10-15 Matsubara frequencies) of p(iu>k) to a 6th- 
order polynomial. The resulting 6 coefficients are then 
used to obtain a (3, 3) Pade approximant using standard 
techniques |2rij |. This approximant is then used for the 
analytic continuation of p by iuok —> lo + ib". Resulting 
real frequency conductivities ct'(uj) are displayed in Fig. 
0Ji versus to /u> c . The typical SSE data are plotted versus 
w/10 and are only shown for L = 20, (3 = 10). 

The results for Hy show a broadened peak as u> — » 0, 
due to inelastic scattering, followed by a second peak 
nicely consistent in height and width with the SSE data. 
The SSE data also displays a high narrow peak at very 
low frequencies, whose position and shape are unstable 
with respect to the choice of the initial image and Max- 
Ent parameters. This is clearly an artifact of the method, 
however its presence is indicative of the tendency to accu- 
mulate the weight at very low frequencies, in qualitative 
agreement with Tiy result. The subsequent fall-off in the 
conductivity at high frequencies is physically consistent, 
but its functional form depends on the Pade approxi- 
mant used. For lo/lo c > 1/2, we expect the analytic 
continuation of the data for Ti.y to become sensitive to 
the order of the approximant used and we therefore in- 
dicate the results in this regime by dotted lines only. We 
note that results at all temperatures yield the same dc 
conductivity a* — 0.45(5)<tq, theoretically predicted Q 



to be universal. Due to the very different scaling pro- 
cedure this result differs from previous numerical result 
(7* = 0.285(20)ctq on the same model in the T = 
limit. It also differs significantly from a theoretical esti- 
mate 0, a* = 1.037(Tq, valid to leading order in e = 3— of. 
Remarkably, our result for the dc conductivity is very 
close to the one obtained in Ref.^for the phase transition 
in the disordered Bose-Hubbard model. Experimental re- 
sults indicate a value close to unity |27| , however it was 
previously observed |7j that long-range Coulomb interac- 
tions, impossible to include in the present study, tend to 
increase a considerably. The same data are shown versus 
lu/T in Fig. 0b. Notably, when using this parametriza- 
tion uj c cancels out and all our data follow the same func- 
tional form. The scaling with uj/T at low frequencies is 
now immediately apparent, with a surprisingly wide low 
uj/T peak. The width of this peak is consistent with the 
data in Fig. [311. Furthermore, on the same lo/T scale the 
continuous time SSE data for Hbh and the results for 
Hv qualitatively agree. 

In summary, we have demonstrated that by doing a 
very careful data analysis it is possible to observe the 
theoretically predicted universal w/T-scaling at the 2D 
superfluid-insulator transition. We have also estimated 
the universal dc conductivity at this transition and found 
that it differs significantly from existing numerical and 
theoretical estimates. 
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